load('storeA.mat','AA');

%Count for end clockwise(1) or anti-clockwise(2)
Count=zeros(2,1);

%Number of experiment
ExpNum=3000;

%Set the starting point
%The starting time is 36.2hr
StartPosTemp=(find(AA(:,2)>36.2));
StartPos=StartPosTemp(1);

TimeDisP=zeros(ExpNum,1);
TimeDisN=zeros(ExpNum,1);

for Single=1:ExpNum
    if mod(Single,100)==0
        display(Single);
    end;
    state=1;
    cycle=0;
    flag=1;
    while flag==1
        dt=AA(StartPos+cycle+1,2)-AA(StartPos+cycle,2);
        prestate=state;
        switch state
            case 1
                Kt2=KCalculation(1,2,AA(StartPos+cycle,1));
                Kt4=KCalculation(1,4,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt2)
                    state=2;
                end;
                if rand(1)<(dt*Kt4)
                    state=0;
                end;
            case 2
                Kt3=KCalculation(2,3,AA(StartPos+cycle,1));
                Kt1=KCalculation(2,1,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt3)
                    state=3;
                end;
                if rand(1)<(dt*Kt1)
                    state=1;
                end;
            case 3
                Kt4=KCalculation(3,4,AA(StartPos+cycle,1));
                Kt2=KCalculation(3,2,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt4)
                    state=4;
                end;
                if rand(1)<(dt*Kt2)
                    state=2;
                end;
            case 0
                Kt1=KCalculation(4,1,AA(StartPos+cycle,1));
                Kt3=KCalculation(4,3,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt1)
                    state=1;
                end;
                if rand(1)<(dt*Kt3)
                    state=-3;
                end;
            case -3
                Kt4=KCalculation(3,4,AA(StartPos+cycle,1));
                Kt2=KCalculation(3,2,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt4)
                    state=0;
                end;
                if rand(1)<(dt*Kt2)
                    state=-2;
                end;
            case -2
                Kt3=KCalculation(2,3,AA(StartPos+cycle,1));
                Kt1=KCalculation(2,1,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt3)
                    state=-3;
                end;
                if rand(1)<(dt*Kt1)
                    state=-1;
                end;
            case -1
                Kt4=KCalculation(3,4,AA(StartPos+cycle,1));
                Kt2=KCalculation(3,2,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt4)
                    state=-4;
                end;
                if rand(1)<(dt*Kt2)
                    state=-2;
                end;
        end;
        
%         if prestate~=state
%             state
%         end;
        
        if abs(state)==4
            if state==4
                Count(1)=Count(1)+1;
                TimeDisP(Single)=AA(StartPos+cycle,2);
            else
                Count(2)=Count(2)+1;
                TimeDisN(Single)=AA(StartPos+cycle,2);
            end;
            flag=0;
        end;
        cycle=cycle+1;
        if cycle>(length(AA)-StartPos-10);
            flag=0;
        end;
    end;
end;

figure; bar(Count);
axis([0 3 0 ExpNum]);
xlabel('clockwise or anti-clockwise');
ylabel('Number of cases');

TimeDisP(find(TimeDisP==0))=[];
TimeDisP=TimeDisP-36.2;
figure; hist(TimeDisP,40);
xlabel('Dwell Time / h');
ylabel('Number of cases');

TimeDisN(find(TimeDisN==0))=[];
TimeDisN=TimeDisN-36.2;
figure; hist(TimeDisN,40);
xlabel('Dwell Time / h');
ylabel('Number of cases');
ylabel('Number of cases');